***
      PROGRAM popbin
***
*
* Evolves a population of binaries using input parameters 
* read from input file binaries.in (M1, M2, P, e, Z, Tmax). 
*
***
      implicit none
*
      INCLUDE 'const_bse.h'
*
      integer i,j,k,jj,nm1,Kstage
      integer kw,kw2,kwx,kwx2,kstar(2)
      integer i1,i2,kdum
*
      real*8 Rsun,CG,sep,Mtotal,tbx1,eccx1,mx3
      real*8 m1,m2,tmax
      real*8 mass0(2),mass(2),z,zpars(20)
      real*8 epoch(2),tms(2),tphys,tphysf,dtp
      real*8 rad(2),lum(2),ospin(2)
      real*8 massc(2),radc(2),menv(2),renv(2)
      real*8 sep0,tb0,tb,ecc0,ecc,aursun,yeardy,yearsc,tol
      PARAMETER(aursun=214.95d0,yeardy=365.25d0,yearsc=3.1557d+07)
      PARAMETER(tol=1.d-07)
      real*8 t1,t2,mx,mx2,tbx,eccx,t3,t
*
************************************************************************
* BSE parameters:
*
* neta is the Reimers mass-loss coefficent (neta*4x10^-13: 0.5 normally). 
* bwind is the binary enhanced mass loss parameter (inactive for single).
* hewind is a helium star mass loss factor (1.0 normally).
* alpha1 is the common-envelope efficiency parameter (1.0).  
* lambda is the binding energy factor for common envelope evolution (0.5).
*
* ceflag > 0 activates spin-energy correction in common-envelope (0). 
* tflag > 0 activates tidal circularisation (1).
* ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0). 
* wdflag > 0 uses modified-Mestel cooling for WDs (0). 
* bhflag > 0 allows velocity kick at BH formation (0). 
* nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1). 
* mxns is the maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1). 
* idum is the random number seed used by the kick routine. 
*
* Next come the parameters that determine the timesteps chosen in each
* evolution phase:
*                 pts1 - MS                  (0.05) 
*                 pts2 - GB, CHeB, AGB, HeGB (0.01)
*                 pts3 - HG, HeMS            (0.02)
* as decimal fractions of the time taken in that phase.
*
* sigma is the dispersion in the Maxwellian for the SN kick speed (190 km/s). 
* beta is wind velocity factor: proportional to vwind**2 (1/8). 
* xi is the wind accretion efficiency factor (1.0). 
* acc2 is the Bondi-Hoyle wind accretion factor (3/2). 
* epsnov is the fraction of accreted matter retained in nova eruption (0.001). 
* eddfac is Eddington limit factor for mass transfer (1.0).
* gamma is the angular momentum factor for mass lost during Roche (-1.0). 
*
      neta = 0.5
      bwind = 1.0
      hewind = 0.0
      alpha1 = 1.0
c      lambda = 1.0 
      ceflag = 3
      tflag = 1
      ifflag = 0 
      wdflag = 1 
      bhflag = 1
      nsflag = 1
      mxns = 2.5
      pts1 = 0.05
      pts2 = 0.01
      pts3 = 0.02
      sigma = 265.0
      beta = 0.125
      xi = 1.0 
      acc2 = 1.5
      epsnov = 0.001
      eddfac = 1.0
      gamma = -1.0
*
* Set the seed for the random number generator. 
*
      idum = 3234
      if(idum.gt.0) idum = -idum
*
* Set the collision matrix.
*
      CALL instar
*
* Open the input file - list of binary initial parameters. 
*
      OPEN(10,file='binaries.in',status='old')
      nm1 = 8.0e6
*
* Open the output files. 
*
c      OPEN(11,file='binaries.out',status='unknown')
c      OPEN(11,file='wr.dat',status='unknown')
c      OPEN(12,file='wd.dat',status='unknown')
      OPEN(13,file='hmxb.dat',status='unknown')
c      OPEN(14,file='bh.dat',status='unknown')
*
      do i = 1,nm1
*

      if(MOD(i,1000).eq.0.0)then 
      print*,i
      endif
* Read in parameters and set coefficients which depend on metallicity. 
*
         READ(10,*)m1,m2,sep,ecc,z
         CALL zcnsts(z,zpars)
*
      Rsun = 6.96e8
      CG= 6.67e-11
      Mtotal = 1.99e30*(m1+m2)
      tb = ((sep*Rsun)**3*39.476/(CG*Mtotal))**0.5/8.64e4
*
         ecc0 = ecc
         tb0 = tb/yeardy
         sep0 = aursun*(tb0*tb0*(mass(1) + mass(2)))**(1.d0/3.d0)
         tb0 = tb
*
* Initialize the binary. 
*
         kstar(1) = 1
         mass0(1) = m1
         mass(1) = m1
         massc(1) = 0.0
         ospin(1) = 0.0
         epoch(1) = 0.0
*
         kstar(2) = 1
         mass0(2) = m2
         mass(2) = m2
         massc(2) = 0.0
         ospin(2) = 0.0
         epoch(2) = 0.0
*
         tphys = 0.0
         tphysf = 13000
         dtp = 0.0
*
* Evolve the binary. 
*
         CALL evolv2(kstar,mass0,mass,rad,lum,massc,radc,
     &               menv,renv,ospin,epoch,tms,
     &               tphys,tphysf,dtp,z,zpars,tb,ecc)
*
* Search the BCM array for the formation of binaries of 
* interest (data on unit 12 if detected) and also output 
* the final state of the binary (unit 11). 

*    Time, 
*    [stellar type, initial mass, current mass, log10(L), log10(r),
*    log10(Teff), core mass, core radius, mass of any convective 
*    envelope, radius of the envelope, epoch, spin, mass loss rate and 
*    ratio of radius to roche lobe radius (repeated for secondary)],
*    period, separation, eccentricity.
*
* In this example we will search for CVs. 
*

         jj = 0
         t1 = -1.0
         t2 = -1.0
         t3 = -1.0
 30      jj = jj + 1
         if(bcm(jj,1).lt.0.0) goto 40
         kw = INT(bcm(jj,2))
         kw2 = INT(bcm(jj,16))
         bcm(jj,30) = yeardy*bcm(jj,30)
*
         i1 = 15
         i2 = 29
         if(kw.gt.kw2)then
            kdum = kw2
            kw2 = kw
            kw = kdum
            i2 = 15
            i1 = 29
         endif 
*for Be WR 



         if(bcm(jj,18).GE.0.6d0.and.bcm(jj,18)
     &   .LE.2.2d1.and.
     &   bcm(jj,30).le.10000.and.bcm(jj,32).ge.0.d0
     &    .and.bcm(jj,32).lt. 1.0d0) then !............


*for Be WD 


*for Be NS

          if(kw2.eq.13)then

            if(kw.le.1.and.bcm(jj,i1).le.1.0)then            
               if(t1.lt.0.0)then
                  t1 = bcm(jj,1)
                  kwx = kw
                  kwx2 = kw2
                  mx = bcm(jj,i1-11)
                  mx2 = bcm(jj,i2-11)
                  tbx = bcm(jj,30)
                  eccx = bcm(jj,32)
c                  WRITE(13,112)m1,m2,tb0,mx,mx2,tbx,eccx,t1
               endif
            end if
              

               if(bcm(jj,i1).ge.1.0)then 
                  if(t2.lt.0.0)then
                   t2 = bcm(jj,1)
                   tbx1 = bcm(jj,30)
                   eccx1 = bcm(jj,32)
                   mx3 = bcm(jj,i1-11)               
                    if (t1.gt.0.0)then
            WRITE(13,112)m1,m2,tb0,mx,mx2,tbx,eccx,t1,t2
                   goto 40
                    endif
                  endif
                endif

            endif

*for Be BH



         
         end if
         goto 30
 40      continue
*
      enddo
*
C 111  FORMAT(f10.1,2i3,3f8.3,1p,e14.6)
 112  FORMAT(2f10.2,3f10.2,f10.2,f10.2,6f10.2)

*
************************************************************************
*
      STOP
      END
***
